#include <assert.h>
#include <mpi.h>
#include <stdlib.h>
#include <iostream>
#include <papi.h>
using namespace std;

typedef unsigned long int INDEX;

double init() { return 0.0;}
double accum(double t, double d) { 
  if(t<d)
    return d;
  else
    return t;
}
double combine(double left, double right) {
  if(left<right)
    return right;
  else
    return left;
}
double scan(double t, double s) { 
  if(s>=t)
    return 1.0;
  else
    return -1.0;
}

int main(int argc, char* argv[]){
  int np;
  int rank;
  int ierr;
  MPI_Status status;
  INDEX INPUT_SIZE;
  INDEX size;

  double myTally;
  double lTally;
  double myTallyt;
  double pTally;
  int stride = 1;
  double* data;
  double* result;
  double* result2;
  long long start1;
  long long end1;
  long long start2;
  long long end2;
  // MPI init
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &np);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  // Load file & init
  if(rank==0){
    INPUT_SIZE=10000;
    data = (double*) new double [INPUT_SIZE];
    result = (double*) new double [INPUT_SIZE];
    result2 = (double*) new double [INPUT_SIZE];
    srand(time(NULL));
    for(INDEX i=0; i<INPUT_SIZE; i++)
      data[i] = ((rand()%((i+1)*INPUT_SIZE))+0.1)/(i+1.0);
  }

  if (rank==0){
    // initialize PAPI library
    unsigned long chkflg = PAPI_VER_CURRENT;
    PAPI_library_init(chkflg);
    if (chkflg != PAPI_VER_CURRENT) {
      cout<<"Error PAPI Library out of date"<<endl;
      exit(1);
    }
  }
  // start timing
  MPI_Barrier(MPI_COMM_WORLD);
  if (rank==0){
    start1 = PAPI_get_real_usec();
  }
  // Broadcast input size
  if (MPI_Bcast(&INPUT_SIZE, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD) != MPI_SUCCESS) {
    cout<<"MPI_Bcast failed"<<endl;
    fflush(stdout);
    MPI_Abort(MPI_COMM_WORLD, -1);
  }

  // Allocate memory for local process
  size = INPUT_SIZE/np;
  double* ldata = (double*) new double [size];
  double* lresult = (double*) new double [size];

  // Distribute data for local process
  ierr = MPI_Scatter(data, size, MPI_DOUBLE,
      ldata, size, MPI_DOUBLE,
      0, MPI_COMM_WORLD);

  // Local process
  myTally = init();
  for(INDEX i=0; i<size; i++){
    myTally = accum(myTally, ldata[i]);
  }

  // Going up
  while(stride < np) {
    //lTally[rank+stride] = myTally[rank];
    //myTally[rank] = combine(myTally[rank], myTally[rank+stride]);
    MPI_Barrier(MPI_COMM_WORLD);
    if (rank%(2*stride) == 0) {
      MPI_Send(&myTally, 1, MPI_DOUBLE, rank+stride, 0, MPI_COMM_WORLD);
      MPI_Recv(&myTallyt, 1, MPI_DOUBLE, rank+stride, 0, MPI_COMM_WORLD, &status);
      myTally = combine(myTally, myTallyt);
    }
    if ((rank-stride)%(2*stride) == 0) {
      MPI_Recv(&lTally, 1, MPI_DOUBLE, rank-stride, 0, MPI_COMM_WORLD, &status);
      MPI_Send(&myTally, 1, MPI_DOUBLE, rank-stride, 0, MPI_COMM_WORLD);
    }
    stride = 2*stride;
  }

  if (rank == 0) {
    myTally = init();
  }

  // Propogate down
  stride = np/2;
  while (stride >= 1) {
    MPI_Barrier(MPI_COMM_WORLD);
    //myTally[rank+stride] = combine(myTally[threadid], lTally[rank+stride]);
    if (rank%(2*stride) == 0) {
      MPI_Send(&myTally, 1, MPI_DOUBLE, rank+stride, 0, MPI_COMM_WORLD);
    }
    if ((rank-stride)%(2*stride) == 0) {
      MPI_Recv(&myTallyt, 1, MPI_DOUBLE, rank-stride, 0, MPI_COMM_WORLD, &status);
      myTally = combine(myTallyt, lTally);
    }
    stride = stride/2;
  }
  MPI_Barrier(MPI_COMM_WORLD);
  pTally = myTally;
  for(INDEX i=0; i<size; i++) {
    pTally=accum(pTally, ldata[i]);
    lresult[i] = scan(pTally, ldata[i]);
  }
  if(rank==0){
    for(INDEX i=0; i<size; i++) {
      result[i] = lresult[i];
    }
    for(int i=1; i<np;i++){
      MPI_Recv(result+i*size, size, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, &status);
    }
  }else{
    MPI_Send(lresult, size, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
  }

  // timing again
  MPI_Barrier(MPI_COMM_WORLD);
  if (rank==0){
    end1 = PAPI_get_real_usec();
  }

  // Clean up
  delete [] ldata;
  delete [] lresult;
  // timing again
  MPI_Barrier(MPI_COMM_WORLD);
  if (rank==0){
    start2 = PAPI_get_real_usec();
    //compare
    double ht=0;
    for(INDEX i=0; i<INPUT_SIZE; i++) {
      if(data[i]>=ht){
        ht=data[i];
	result2[i]=1.0;
	//assert(result[i]==result2[i]);
      }else{
	result2[i]=-1.0;
	//assert(result[i]==result2[i]);
      }
    }
    end2 = PAPI_get_real_usec();
    cout<<np<<"p "<<end1-start1<<" s "<<end2-start2<<" size "<<INPUT_SIZE<<" # "<<np<<endl;

  }
  MPI_Finalize();
  return 0;
}
